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state Potts model are performed. The phase diagram and the physical properties at 
the phase transitions are studied using finite-size scaling techniques. Evidences are 
given for the existence of a tricritical point dividing the phase diagram into a regime 
where the transitions remain of first order and a second regime where the transitions 
are softened to continuous ones by the influence of disorder. In the former regime, 
the nature of the transition is essentially clarified through an analysis of the energy 
probability distribution. In the latter regime critical exponents are estimated. Rare 
and typical events are identified and their role is qualitatively discussed in both 
regimes. 
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1 Introduction 



The influence of disorder is of great interest in physics, since pure systems are rare 
in nature. It has been known for more than thirty years that the universahty class 
associated with a continuous phase transition can be changed by the presence of quenched 
impurities [Tj. According to the Harris criterion j^, uncorrelated randomness coupled 
to the energy density can only affect the critical behaviour of a system if the critical 
exponent a describing the divergence of the speciflc heat in the pure system is positive. 
This has been established in the case of the g-state Potts model in dimension D = 2 
for example. For 2 < q < 4, the pure system undergoes a continuous transition with 
a positive critical exponent a. As predicted by the Harris criterion, new universality 
classes have been observed both perturbatively and numerically 3 (for a review, see 
Ref. [1]). The special case q = 2, the Ising model, is particularly interesting since in 
the pure system, the specific heat displays a logarithmic divergence (a = 0) making the 
Harris criterion inconclusive Based on perturbative and numerical studies, it is now 
generally believed that the critical behaviour remains unchanged apart from logarithmic 
corrections when introducing randomness in the system P|. In three dimensions (3D), 
the disordered Ising model was subject of really extensive studies (see, e.g., Ref. [7] for 
an exhaustive list of references). 

Less attention has been paid to first-order phase transitions. It is known that ran- 
domness coupled to the energy density softens any temperature-driven first-order phase 
transition [Hj. Moreover, it has been rigorously proved ^ that in dimension D < 2 
an infinitesimal amount of disorder is sufficient to turn any first- order transition into 
a continuous one. The first observation of such a change of the order of the transi- 
tion was made in the 2D 8-state Potts model \W\ where a new universality class was 
identified jlll I12j . For higher dimensions, the first-order nature of the transition may 
persist up to a finite amount of disorder. A tricritical point at finite disorder between 
two regimes of respectively first-order and continuous transitions is expected jlll I13j . 
The existence of such a tricritical point for the site-diluted 3D 3-state Potts model could 
only be suspected by simulations because the pure model already undergoes a very weak 
first-order phase transition On the other hand, the first-order phase transition of 
the pure 5-state Potts model is very strong and would hence make it rather diflicult to 
study the role of disorder. As a consequence, we have turned our attention to the 3D 
4-state Potts model and have shown that there exists a second-order transition regime 
for this model JH]- Our choice of bond dilution is motivated by the fact that for this 
model only high-temperature expansions results are available up to now which to our 
knowledge cannot be done for site-dilution or are at least more difficult |16j . 

In Sect.ISJwe define the model and the observables, and remind the reader of how these 
quantities behave at first- and second-order phase transitions. Section |31 is devoted to 
the numerical procedure, first the description of and then the comparison between the 
algorithms which are used at low and high impurity concentrations, followed by a first 
discussion of the qualitative properties of the disorder average. A short characterisation 
of the nature of the phase transition - at a qualitative level - is reported in Sect. 01 The 
motivation of this section is to first convince ourselves that the transition does indeed 
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undergo a qualitative change when the strength of disorder is varied. Then, we describe 
how the phase diagram is obtained and concentrate on the first-order regime in Sect.|21 
In Sect-ini we discuss the critical behaviour in the second-order induced regime. Finally, 
the main features of the paper are summarised in Sect. [Tj 



2 Model and observables 

We study the disordered 4-state Potts model on a cubic lattice A. The model is defined 
by the Hamiltonian 

n[a,J] = -Y,J,,6„,^,„^, (1) 

where the spins ctj, located on the vertices i of the lattice A, are allowed to take one 
of the g = 4 values = The boundaries are chosen periodic in the three 

space directions. The notation 7^[o", J] specifies that the Hamiltonian is defined for any 
configuration of spins and of couplings. The sum runs over the couples of nearest- 
neighbouring sites and the exchange couplings Jij are independent quenched, random 
variables, distributed according to the normalised binary distribution (J > 0) 

P[M = X{\pKJ'.i -J) + {i-pW^j)]. (2) 

The pure system (at p = 1) undergoes a strong first-order phase transition with a 
correlation length ^ ~ 3 lattice units at the inverse transition temperature [3c J = 
0.62863(2) J7j (we keep the conventional notation j3 = {ksT)^^, since in the con- 
text there is no risk of confusion with the critical exponent of the magnetisation). As 
far as we know, no more information has been made available on this model. We do 
not expect any phase transition for bond concentration p smaller than the percolation 
threshold pc = 0.248 812 6(5) ^H] since the absence of a percolating cluster makes the 
appearance of long-range order impossible. 

In the following, we are thus dealing with quenched dilution. The averaging pre- 
scription is such that the physical quantities of interest in the diluted system (say an 
observable Q) are obtained after averaging first a given sample [J] over the Boltzmann 
distribution, {Q[j])/3, and then over the random distribution of the couplings denoted 

by {Q[j])i3, since there is no thermal relaxation of the degrees of freedom associated to 
quenched disorder: 

i> The thermodynamic average of an observable Q at inverse temperature /? and for 
a given disorder realization [j] is denoted 

(Q[j])/3 = (%](/3))-^/pHg[.,j]e-'3«['^'^i 
1 
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where Nmcs is the number of Monte Carlo iterations (Monte Carlo steps) during 
the "production" part after the system has been thermalised. Here we use the 
following notation: Q[c7,j] is the value of Q for a given spin configuration [a] and a 
given disorder realization [J], is a value obtained by Monte Carlo simulation 

at inverse temperature /3. Time to time, we will have to specify a particular disorder 
realization, say #n, and the value of the observable Q for this very sample will be 
denoted as 

> The average over randomness is then performed, 



{QlJ])l3 



V[J]{Q[j])pP[J] 



N{J} 



[J] 

dQQPpiQ), 



(4) 



where N{J} is the number of independent samples. The probability P/3{Q) is 
determined empirically from the discrete set of values of {Q[j])i3. This disorder 
average is simply denoted as Q{P) for short, i.e. Q{f3) = {Q[j])f3- 

For a specific disorder realization [J], the magnetisation per spin fn^^^j] = L^^M^^ j] 
of the spin configuration [a] is defined from the fraction of spins, P[a,j], that are in the 
majority orientation. 



max 

O"0 



QP[a,J] - 1 

q-1 



(5) 



The order parameter of the diluted system is thus denoted m(/?) = (m[j])/3. Thermal 
and disorder moments (mj^j)^ and respectively, are also quantities of interest. 

The magnetic susceptibility X[J]{P) and the specific heat C[j](/3) of a sample are defined 
using the fiuctuation-dissipation theorem, i.e. 



X[j]i(3) 
C[j]{f3)/kB 



(6) 
(7) 



where 

e[a,J]=L' E[^j^=L' ^Jij^(T„ay (8) 

is the negative energy density since -E[o-,j] = —'T~i[a, J]. Binder cumulants ^H] take their 
usual definition, for example 



f^m,,„(/3) = l 



3(m2^])/3 



2 • 



(9) 
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Derivatives with respect to the exchange coupUng are computed through 

L YpMm^A)p= ^^.p^ -(^W)/^- (10) 
All these quantities are then averaged over disorder, yielding C{(3), UmiP), and 



At a second-order transition, these quantities are expected to exhibit singularities 
described in terms of power laws from the deviation to the critical point. These power 
laws define the critical exponents. In the following, the properties will be investigated 
using finite-size scaling analyses, i.e., according to the following size dependence at the 
critical temperature, 

m(/3„L-i) ~ B.L-"/'', (11) 
xiPcL-') ~ ^,L^/^ (12) 
C(/3e,L-i) ~ (13) 



dp 



Nn,cL^'r (14) 



At a first-order transition, the order parameter has a discontinuity at the transition 
temperature, suggesting that j3/v formally becomes zero. Heuristic (and for pure q- 
state Potts models with sufficiently large q even rigorous) arguments also suggest that 
a/u, and 1/v should then coincide with the space dimension D \2G., restoring the 
ordinary extensivity of the system in Eqs. H12|) - (|14() . 

When the transition temperature is not known exactly, the problem of the value of 
the inverse critical temperature f3c in the expressions above can be a source of further 
difficulties. Usually, one follows a flow of flnite-size estimates given by the location 
of the maximum of a diverging quantity (for example the susceptibility, Xmax(-^~^) = 
max^[)^(/3, L~^)]). From the scaling assumption, supposed to apply at the random fixed 
point, 

xW,L-^) = L^'-f^{L^'-t), (15) 
with t = |/3c — the inverse temperature /3max where the maximum of x occurs, 

scales according to 

/3max~/3c + a^"'/". (17) 

Notice that the scaling function f^ix) takes its maximum value = Xmax^ '^^'^ 

X = a. 

At that very temperature where the finite-size susceptibility has its maximum, we 
then have similar power law expressions, 

m(/3„,ax,^"') ~ /™(a)L-^/^ (18) 
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~ (19) 
~ /c(a)L"/^ (20) 

~ /m,n(a)Li/^ (21) 

These equations are similar to Eqs. - (|14() where the amplitudes take the values 
Be = /m(0), Te = /x(0), A, = /c(0), and iV, = /^,„(0). 

From this discussion, we are led to give a more precise definition of Xmaxi^^^)- 
reasonable alternative definition of this quantity could be the disorder average of the 
individual maxima corresponding to the different samples. Each of them has its own 
susceptibility curve, xiJjiP^ which displays a maximum X[j],max(-^~^) at a given 
value of the inverse temperature /3™p(L^^). These values X[j],max(-^~"^) then be 
averaged, but this is in general different from the definition that we gave for the average 
over randomness in Eq. (jlj. Here we keep as a physical quantity the expectation value 
x(/3, L~^) which is then plotted against /3, and f3ma.x{L^^) in Eq. (fT^ is the temperature 
where the disorder averaged susceptibility displays its maximum which is thus identified 
with Xmax(-^~^)- In the following, this is the physical content that we understand when 

discussing Xmax(-^^""^)- 

3 Numerical procedures 

We conducted a long-term and extensive study of the bond-diluted 3D 4-state Potts 
model, and it is the purpose of this paper to report results for moderately large system 
sizes in the first-order regime, and an extended analysis based on really large-scale com- 
putationsin the second-order regime. Cross-over effects between different regimes are 
also discussed. The simulations were performed on the significant scale of several years. 
A strict organisation was thus required, and we proceeded as follows: as an output of 
the runs, all the data were stored in a binary format. For each sample (with a given 
disorder realization and lattice size) and each simulated temperature, the time series 
of the energy and magnetisation were stored. A code was written in order to extract 
from all the available files the histogram reweightings of thermodynamic quantities of 
interest, entering as an input the chosen dilution, lattice size, temperature, .... It is also 
possible to adjust the number of thermalisation iterations, the length of the production 
runs where the thermodynamic averages are performed, the number of samples for the 
disorder average, or to pick a specific disorder realization, and so on. In some sense, the 
time series correspond to the simulation of the system, and we can then measure physical 
quantities on it, and virtually produce as many results as we want. Of course, this is 
not what we intend to do in the following, we rather shall try to concentrate only on the 
most important results. 



C(/3max)^ ) 



df3 
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Figure 1: Evolution of the susceptibility as the size of the system increases (up to i = 10) in the 
two different regimes: pure system p = 1.0 on the left plot and high dilution p = 0.48 
on the right plot. 

3.1 Choice of update algorithms 

We studied this system by really large-scale Monte Carlo simulations. A preliminary 
study, needed in order to schedule such large-scale Monte Carlo simulations, showed 
that the transitions at small and high concentrations of non- vanishing bonds p were (as 
expected) qualitatively different: 

i> Close to the pure system, p ~ 1, the susceptibility peaks develop as the size 
increases to become quite sharp (see Fig. ^ , in agreement with what is expected 
at a first-order phase transition 21 . 

> At larger dilutions (small values of p) on the other hand, the peaks are softened 
and are compatible, at least at first sight, with a second-order phase transition 
(see Fig.lIJ). 

As will be demonstrated below, the tricritical dilution dividing these two regimes is 
roughly located at ptcp ~ 0.68 — 0.84. In the regime of randomness-induced continuous 
transitions (or weak first-order transitions, that is at low non-zero bond concentration 
p), the Swendsen-Wang cluster algorithm ^2 preferred in order to reduce the crit- 
ical slowing-down. As already pointed out by Ballesteros et al. |23| I14j. a typical spin 
configuration at low bond concentration is composed of disconnected clusters for most 
of the disorder realisations. It is thus safer to use the Swendsen-Wang algorithm, for 
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which the whole lattice is swept at each Monte Carlo iteration, instead of a single-cluster 
Wolff update procedure. In the strong first-order regime (high bond concentration p), 
the multi-bondic algorithm a multi-canonical version of the Swendsen-Wang algo- 
rithm, was chosen in order to enhance tunnellings between the phases in coexistence at 
the transition temperature. The Swendsen-Wang algorithm, being less time-consuming, 
was nevertheless preferred even in this regime of long thermal relaxation as long as at 
least ten tunnelling events between the ordered and disordered phases could be observed. 
As the first-order regime is approached, more and more sweeps are needed to fulfil this 
condition. We had to use up to 

- 200 000 Monte Carlo steps (MCS) at p = 0.76 for L = 16 for example, 
while in the second-order regime for much larger systems, we needed 

- 100 000 MCS at p = 0.68 for L = 50, 

- 30 000 MCS atp = 0.56 for L = 96, and 

- 15 000 MCS atp = 0.44 for L = 128. 

This is the essential reason for the size limitation in the first-order regime^. A com- 
parison between the two algorithms is illustrated in the case of the pure system for a 
moderate size {L = 6) in Fig. [2^. The insert shows a zoom of the peak in the suscepti- 
bility and reveals as expected that in this first-order regime, the multi-bondic algorithm 
provides a better description of the maximum which, since being higher is probably 
closer to the truth. 

In both regimes, the procedure of histogram reweighting enables us to extrapolate 
thermodynamic quantities to neighbouring temperatures. It leads to a better estimate 
of the transition temperature and of the maximum of the susceptibility, refining the 
finite-size estimate at each new size considered, since the maximum is progressively 
reached (Fig. I^Ja and[2t). The reweighting has to be done for each sample, then the 
average is obtained as in Eq. Q). For a particular sample, the probability to measure 
at a given inverse temperature /3 a microstate [a] with total magnetisation Mj^. j] = M 
and energy Ej^^jj = E, is PpiM^E) = {Z^j^{0))-^Q.{M,E) e^^ where n{M,E) is the 
degeneracy of the macrostate. Note that we defined E as minus the energy in order to 
deal with a positive quantity. We thus get at a different inverse temperature /?' 

P^,{M,E) = {ZyjP)/Zyj0))Pp{M,E) e(/^'-^)^, (22) 

where the prefactor (Z[j](/3)/Z[j](/3')) only depends on the two temperatures. For any 
quantity Q depending only on M\^„^j^ and -E[cr,j] the thermal average at the new point (3' 
hence follows from 

,^ , _ EM,EQiM,E)Pp{M,E)e(^'-^)^ 

A rough estimate of the time needed by a single simulation is given by X (#MCS) X 1/is for one 
sample and one temperature. 
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Figure 2: a) Comparison between canonical Swendsen-Wang and multi-bondic algorithms 
for a pure system (j> — 1.0) of size L — 6 (histogram reweightings 
produced from simulations at inverse temperatures /3J — 0.605 to 0.655 
are superimposed). The insert shows a zoom of the peak location, 
b) and c) Histogram reweighting of the average susceptibility xiP) in a disordered 
system with p — 0.56 at different sizes L — 25 and 64. The maximum is progressively 
obtained after a few iterations (the next simulation is performed at the temperature 
of the maximum of the histogram reweighting of the current simulation). 



9 



It is well known that the quality of the reweighting strongly depends on the number of 
Monte Carlo iterations, the larger this number the better the sampling of the configura- 
tion space and thus of the tails of Pf^. Here we have to face up to the disorder average 
also and a compromise between a good disorder statistics and a large temperature scale 
for the reweighting of individual samples has to be found, but we are mainly interested 
in the close neighbourhood of the susceptibility maximum, i.e., in a small temperature 
window. 



3.2 Equilibration of the samples and thermal averages 

Before any measurement, each sample has to be in thermal equilibrium at the simulation 
temperature. Starting from an arbitrary initial configuration of spins, during the initial 
steps of the simulation process, the system explores configurations which are still strongly 
correlated to the starting configuration. The typical time scale over which this "memory 
effect" takes place is measured by the autocorrelation time. The integrated energy 
autocorrelation time (one can define more generally an autocorrelation time for any 
quantity) is given by 

I — U J — i 

where a'^ = {efj])f3 — is the variance, e-jjj is the value of the energy density at 

iteration j for the realization [J], and / is a cutoff (as defined, e.g., by Sokal |25j ) 
introduced in order to avoid to run a double sum up to A'^mcS) which would render the 
estimates very noisy. 

It is worth giving a definition of the errors as computed in this work. There are two 
different contributions. Assuming the different realizations of disorder as completely 
independent, one has an error due to randomness on any physical quantity Q, defined 
according to 

To the thermal average for each sample is also attached an error which depends on the 
autocorrelation time t^, such that the total error on a physical quantity is here defined 
as: 

A««=(^A?,„0 + _^j . (26) 

For each disorder realization, the preliminary configurations have to be discarded and 
one usually considers that after 20 times the autocorrelation time r^, thermal equilibrium 
is reached. The measurement process can then start and the thermal average of the 
physical quantities is considered, in the case of a single sample, as satisfying when 
measurements were done during typically 10^ x r*^. For a quantity Q, a satisfactory 
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p=0.44 




5000 10000 
iteration # 



15000 



Figure 3: Monte Carlo data of the magnetisation for the disorder reahzation that gave the largest 
value of X[j] (/3max) for p ~ 0.44 at lattice size L ~ 128, p ~ 0.56 {L ~ 96) and p = 0.68 
(L = 50). 




Figure 4: Susceptibility for sample #1 for p = 0.44 and L = 128. The different curves show the 
result of histogram reweighting of simulations close to the maximum location after 
5 000, 10 000, 12 000, and 15 000 MCS. We can safely consider that the value at the 
temperature of the maximum of the average susceptibility (vertical dashed line) is 
reliable after 15 000 MCS. 
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relative error of the order of 



■therm. Q 



2tQ 



(27) 




^ Af{J} X A^MCS 



~ 10 



indeed requires typically N{J} x A^mcs — W'^t^. Since we also need a large number 
of disorder realizations in order to minimise 

^?dm<3' typically N{J} ~ 10^ - 10^ each 
sample requires a "production" process during A^mcs - (10° - 10'^)tQ. In this paper, we 
choose to work at the upper limit with A'^mcs > 10^ r*^ (since there is a single dynamics 
in the algorithm, the time scale is usually measured through the energy autocorrelation 
time) and N{J} > 10^. 

Examples of times series of the magnetisation are shown in Fig.l^lfor particular samples 
(those which contribute the most to the average susceptibility) at the three largest sizes 
studied^ at dilutions p = 0.44, 0.56, and 0.64 in the second-order regime. The simulation 
temperature is extremely close to the transition temperature and tunnelling between 
ordered and disordered phases guarantees a reliable thermal average. 

Another test of thermal equilibration is given by the influence of the number of MCS 
which are taken into account in the evaluation of thermal averages. An example is shown 
in Fig. 0] where the histogram reweightings of the susceptibility, as obtained with different 
MCS #'s, are shown for a typical sample (the first sample, # 1, is in fact supposed to be 
typical). Although quite different far from the simulation temperature (which is close 
to the maximum) the different curves are in a satisfying agreement at the temperature 
/3max of the maximum of the average susceptibility, shown by a vertical dashed line. The 
criterion 7^ MCS > 250 x t'^ is safely satisfied for the larger number of iterations^. 

3.3 Properties of disorder averages 

For different samples, corresponding to distinct disorder realizations, the susceptibility 
X[.j]iP) at thermal equilibrium may have very different values (see Fig. |21 where the 
running average over the samples is also shown and remains stable after a few hundreds 
of realizations). 

We paid attention to average the data over a sufficiently large number of disorder 
realizations (typically 2000 to 5000) to ensure reliable estimates of non-self-averaging 
quantities [2^]. Averaging over a too small number of random configurations leads to 
typical (i.e. most probable) values instead of average ones. Indeed, as can be seen in 
Fig. El the probability distribution of X[j] (plotted at the temperature /3max where the 
average susceptibility is maximum) presents a long tail of rare events with large values 
of the susceptibility. These samples have a large contribution to the average, shifted far 
from the most probable value. The larger the value of p, the longer the tail. Scanning 
the regime close to the first-order transition thus requires large numbers of samples to 

^The main illustrations are shown in the worst cases, i.e., for the largest systems at each dilution. 
^We also note that something happened between 5000 and 10000 MCS, since the shape of x#i at high 

temperatures becomes unphysical. This is an illustration of the finite window of confidence of the 

histogram reweighting procedure. 
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Fi gure 5: Different values of X[J](/5max) (over the saraples) of susceptibility at p = 0.56, L = 
96 (the simulation is performed at the temperature of the maximum of the average 
susceptibility). The running average over the samples is shown by the solid line. 

explore efficiently the configuration space, so tlie simulations were limited to L = 50 at 
p = 0.68 while we made the calculation up to L = 128 at p = 0.44. In the example of 
Fig. 121 the thermodynamic quantities have been averaged over 3500 disorder realizations 
for p = 0.44 at lattice size L = 128, 2048 for p = 0.56, L = 96, and 5000 disorder 
realizations for p = 0.68, L = 50. 

Self-averaging properties are quantified through the normalised squared width, for 
example in the case of the susceptibility, = c7^(L)/x^, where cr^ = ~x'^- For 
a self-averaging quantity, say Q, the probability distribution, albeit not truly Gaus- 
sian, may be considered so in first approximation close to the peak, and P/siQ) — 

(27r(TQ)~"'^/^e~^'^~^'^^'* ^^'^Q evolves towards a sharp peak in the thermodynamic limit, 

PfiiQ) -^L~*oo S{Q — {Q))- The probability of the average event (Q) goes to 1 and 
the normalised squared width evolves towards zero in the thermodynamic limit while it 
keeps a finite value for a non-self-averaging quantity, as shown in the case of the suscep- 
tibility in Fig. [3 The observation of a longer tail in the probability distribution of the 
X values when p increases is expressed in Fig. [7|by the fact that x becomes less and less 
self-averaging when p increases. 

In contradistinction to the magnetic susceptibility, the energy seems to be weakly self- 
averaging in the range of lattice sizes that we studied as seen in Fig. |H1 The associated 
exponent depends on the concentration of bonds p. This concentration dependence may 
be effective and due to corrections generated by other fixed points (see below). 

In Table Q the influence of the number of MCS is shown for typical samples, but 
also for the average susceptibility. Although the variations for a given sample and from 
sample to sample are important, the average seems stable with our choice of number of 
iterations (the largest), and also the autocorrelation time (for the average) is stable. 
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Figure 6: Probability distribution of the susceptibility X[./] (/^max) for the bond concentrations 
p = 0.44, 0.56, and 0.68 for the largest lattice size in each case. The full curve 
represents the integrated distribution. At each dilution, a full vertical line shows the 
location of the average susceptibility, a dashed line shows the median and a dotted 
line shows the average over the events which are smaller than the median. 



Table 1: Evolution of the susceptibility with the number of Monte Carlo sweeps per spin for 
different samples, X[J] ^-nd the average value (with 2048 samples) ai p — 0.56, L = 96. 
The data are given at the maximum location of the average susceptibility, /3max- The 
last column gives the number of independent measurements per sample. 



A'mcs 


X#i 


X#2 


X#3 


X#4 


X#5 


Xmax 


r^(/3max) 


meas. / sample 


5000 


994 


404 


611 


682 


1803 


617(8) 


95.1 


~ 50 


10000 


952 


390 


698 


614 


1574 


634(8) 


107.4 


~ 90 


15000 


1010 


356 


680 


819 


1398 


638(8) 


111.7 


~ 130 


20000 


939 


351 


689 


851 


1320 


641(7) 


114.0 


~ 175 


25000 


911 


327 


675 


848 


1308 


643(8) 


115.3 


~ 200 


30000 


934 


327 


733 


837 


1297 


643(8) 


116.9 


~ 250 
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Figure 7: Normalised squared width of the susceptibihty, i?^, plotted against the inverse lattice 
size for the three dilutions p = 0.44, 0.56, and 0.68. The solid lines are polynomial 
fits used as guides for the eyes. Note that x is apparently less and less self-averaging 
as p increases. 

In Fig. El a full vertical line points out the location of the average susceptibility Xmax- 
In order to give a comparison, the median value Xmed) defined as the value of X[.j] where 
the integrated probability takes the value 50%, is shown as the dashed line. The more 
it differs from the average, the more asymmetric is the probability distribution. This is 
more pronounced when p increases. We also notice that the maximum of the probability 
distribution (the typical samples) corresponds to smaller susceptibilities. For a given 
number of disorder realizations, this peak is better described than the tail at larger 
susceptibilities, so we also define (shown as dotted lines) an average over the samples 
smaller than the median susceptibility, that we denote X50%) 



where the factor 2 normalises the truncated distribution. In the particular case of the 
probability distributions observed here, i.e. with a sharp initial increase, a peak located 
at small events and a long tail at large values of the variable^, this definition empirically 
gives a sensitive measure of the typical or most probable value. We shall refer to this 
quantity when typical behaviour will be concerned. 



*This shape of probability distribution is very different than in the case of the 3D dilute Ising model |27|. 




(28) 
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Figure 8: Normalised squared width of the energy, Re plotted on a log- log scale against the 
inverse lattice size for the three dilutions p = 0.44, 0.56, and 0.68. Power-law fits 
have been performed and corresponding exponents printed by the curve. The insert 
presents the same data plotted on a linear scale. 

4 Qualitative description of the transition 

Before performing a quantitative analysis of the transition, it is interesting to study in 
some detail why the probability distributions have significantly different shapes when 
p varies, and which type of sample can be considered as a typical one, or which one 
corresponds to a rare event with quite a large or very small susceptibility. Here we shall 
focus on the second-order regime and in particular on p = 0.44 for the largest simulated 
size, L = 128, for which the probability distribution of X[j] at /3max can be inspected in 
Fig. El 

Each sample displays its own maximum and due to the fluctuations over disorder, the 
temperature /?[j]^max where it occurs varies from sample to sample. In Tabled we quote 
for a few rare and typical samples the values of X[J],nax /?[J],max) the maximum of the 
sample susceptibility and the corresponding inverse temperature (see also Figs. inito[TT|). 
The relative variations of these numbers with respect to their average values at /3max : 

^X[J]/Xmax = [X[J],max " Xmax]/Xmax, and A/3[j]//3max = [/3[J],max " /?max]//?max are also 

collected in Table|2j It turns out that rare events with large susceptibility do also display 
a very small shift of the temperature /3[j],max with respect to the average. Other events 
have a smaller susceptibility at /3max both because their maximum X[j],max is smaller but 
also because of a larger shift of the temperature /3[j],max where this maximum occurs. 
A few examples of rare events corresponding to large values of X[,j] are shown in Fig. [HJ 
Rare events corresponding to small values of X[j] are presented in Fig. IllL They have 
a very small contribution to the phase transition so in the following, we will refer only 
to events with large values of X[j] when mentioning rare events. In Fig. 1101 the same 
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Table 2: Relative variations of the peak height Ax[j]/xmax and peak location A/3[j] ,„ax//3max 
for a few samples, chosen among the rare and the typical samples at p = 0.44, L — 128. 
For reference, the values of the average are given by /3max>/ — 1-4820, Xmax — 1450. 
The asterisks (*) mark those samples that are discussed in detail in Figs. I12I-I14I 



type 


sample # 


X[Jl,max 


(^\J],nia,xJ 


^X[Jl/Xmax 


^P\J] ,max//5max 




0035 (*) 


5253 


1.4823 


262.3 


% 


0.02 % 


rare 


0438 


3862 


1.4822 


166.3 


% 


0.013 % 


(large x) 


1135 


3825 


1.4821 


163.8 


% 


0.007 % 




3302 


4314 


1.4823 


197.5 


% 


0.02 % 




0006 


1550 


1.4831 


6.9 


% 


0.07 % 


typical 


0008 (*) 


2792 


1.4810 


92.5 


% 


-0.07 % 


(around peak) 


0021 


1473 


1.4819 


1.6 


% 


-0.007 % 




0039 


2345 


1.4817 


61.7 


% 


-0.02 % 




0373 


946 


1.4852 


-34.7 


% 


0.22 % 


rare 


1492 (*) 


286 


1.4830 


-80.3 


% 


0.07 % 


(small x) 


1967 


1063 


1.4847 


-26.7 


% 


0.2 % 




2294 


769 


1.4853 


-46.9 


% 


0.2 % 



is done for typical events, i.e., those for which the values of X[J] in the peak of 
the distribution. The scales of both axis are the same in the three figures in order to 
facilitate the comparison. 

In Figs. 1121141 we can follow the fluctuations of the magnetisation during the ther- 
malisation process (after equilibration) for three different samples. Configuration ^^35 
(Fig. I12j) corresponds to a rare event, with the definition given above, while the sample 
#8 (Fig. [T^ is a typical one. The last sample, #1492 (Fig. [Tl]) . is an example of a 
realization of disorder which leads to a very small susceptibility peak. These figures also 
present the magnetisation and energy probability distributions. The rare event (Fig. 112(1 
displays a double-peak structure in the probability distributions (only a shoulder is vis- 
ible in -f/3jnax(^[J]))> presumably a remnant of the first-order type transition of the pure 
system. In the average behaviour, it seems that at small values of p, these types of 
samples are "lost" in the large majority of typical samples which have a "second-order 
type" of probability distribution. This observation is corroborated by similar "signals" 
in Figs. El to ^2 concerning the shape of the susceptibility (narrow peak for rare events 
with large susceptibilities and broader for others), of the order parameter (sharp increase 
with P at the transition for rare events, and smoother variation for the typical samples), 
or of the Binder cumulant (deep well at the transition in the case of rare events and less 
pronounced wells for typical ones). 

We may thus argue that a possible mechanism which keeps the pure model's first- 
order character of the transition at larger values of p is connected to the occurrence of 
a larger proportion of samples with the "first-order type", i.e., a very big susceptibility 
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1.480 1.485 




1.480 1.485 




1.480 



1.485 



Figure 9: Examples of rare events for p = 0.44 and L — 128. with large values of X[j]- The 
thick lines show the averages over all realizations. 



signal at /3max- In Fig. El the quite long tail of large susceptibilities in the susceptibility 
distribution confirms this assumption, for p = 0.84 [L = 20), i.e., closer to, or probably 
inside the first-order regime. Also the double-peak structure of the energy distribution at 
this dilution (see Fig. I16() is compatible with a first-order like transition for the average 
behaviour (of course one would have to study the evolution of the energy barrier as the 
size increases, but this makes no sense for a specific disorder realization for which the 
notion of thermodynamic limit is meaningless). The possible interpretation is that the 
rare events of higher susceptibilities when p becomes larger are more comparable to a 
system displaying a first-order transition. This would explain that the susceptibility 
peak is narrower (and thus does coincide with the temperature of the maximum of the 
average only in very rare cases). 



5 Phase diagram and strength of the transition 
5.1 Transition line 

We can now come back to the preliminary phases of this work. The transition temper- 
ature was determined for 19 values of the bond concentration ranging from p = 0.28 to 
p = 1.00 (pure system). We defined an effective inverse transition temperature l3c{L,p) 
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Figure 10: Examples of typical events for the same parameters as in Fig. El The thick lines 
show the averages over all realizations. 



at a given lattice size L as the location of the maximum of the average magnetic suscep- 
tibility X (see Fig. ITT]) . Any diverging quantity could equally have been chosen but it 
turned out that the specific heat was displaying larger statistical errors than the mag- 
netic susceptibility. Moreover, the stability of the random fixed point implies a slowly 
varying specific heat with a critical exponent a < 0.^ 

For each p and L, several Monte Carlo simulations were necessary to get a reasonable 
estimate of PciL,p)- As mentioned before, histogram reweighting was used to refine the 
determination. The procedure was applied up to lattice sizes L = 16. The resulting 
phase diagram for two different lattice sizes is plotted in Fig. 1181 The data appear to be 
in a remarkable accordance. 

The numerical data presented in Fig. ^1 are furthermore in agreement with the mean- 
field prediction Tc{p) = pTc{p = 1) for large bond concentration, close to the pure 
system {p — 1). At smaller concentration p, the topological properties of the bond 
configuration become important and the mean-field prediction fails to reproduce the 
observed behaviour. The effective-medium approximation introduced in this context in 
the eighties by Turban j2H] reproduces quite accurately the numerical data. Limiting the 



^We expect a stable randomness fixed point at large enough dilutions, where the exponent a should be 
negative hence the singular contribution to the specific heat would not be diverging. 
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Figure 12: Time scries of the magnetisation and the energy density and corresponding probabil- 
ity distributions for a rare event (#35) with large X[j] {p = 0.44, L = 128, simulation 
at inverse temperature (3J = 1.48218). 
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Figure 13: Time series of the magnetisation and the energy density and corresponding probabil- 
ity distributions for a typical event (#8) with X[J\ in the vicinity of the most probable 
value {p = 0.44, L = 128, simulation at inverse temperature jSJ = 1.48218). 
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Figure 14: Time series of the magnetisation and the energy density and corresponding probabil- 
ity distributions for a rare event (#1492) with very small X[j], which looks similar to 
typical events {p = 0.44, L = 128, simulation at inverse temperature /3J = 1.48218). 
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Figure 15: Probability distribution of the susceptibility for a system of size L — 20 at p = 0.84 
and PJ = 0.74704, in the seemingly first-order regime. The simulation was performed 
with the multibondic algorithm. 
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Figure 16: Probability distributions P and P[j] of the energy e for the average behaviour and 
for a rare event (large susceptibility), respectively, at p = 0.84 (L = 13). The 
double-peak structure suggests a behaviour for this specific sample which is similar 
to the one observed at a first-order transition. The simulation is performed at inverse 
temperature (3 J = 0.746 356. 
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Figure 17: Average susceptibility and its histogram reweighting for systems of sizes 10^ and 16^ 
for dilutions (from left to right) p = 0.32, 0.40, 0.48, 0.56, 0.64, 0.72, 0.80, 0.88, and 
0.96. 



approximation to a single bond, the following estimate for the transition temperature is 
obtained: 



/3e(p) = J-Mn 



(29) 



{P - Pc) 

where /3P"''°J = 0.62863(2) for the pure system. This expression is exact (as exact as 
it might be with numerical factors introduced) in the limits of the pure system (p = 1) 
and the percolation threshold {pc = 0.248812 6(5)). 



5.2 Order of the transition 

Distinguishing a weak first-order phase transition from a continuous one is a very difficult 
task. The autocorrelation time of the energy at the transition temperature may be 
useful, since it displays a behaviour which depends on the order of the transition. When 
using a canonical Monte Carlo simulation for the study of a first-order transition, the 
time-scale of the dynamics is dominated by the tunnelling events between the ordered 
and disordered phases in coexistence at the transition temperature. Such a tunnelling 
event implies the creation and the growth of an interface whose energy cost behaves as 
PAF = 2ao.d.L^~^ where ao.d. is the reduced interface tension. As a consequence, the 
autocorrelation time grows exponentially as 

r^(L) -e^'^-d-^''"'. (30) 

For a continuous phase transition, this interface tension vanishes and the autocorrelation 
time scales as a power-law of the lattice size, 

r%L) ~ L^ (31) 
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Figure 18: Transition temperatures kBTc{p)/J with respect to the bond concentration p for two 
lattice sizes L — 10 and L = 16. Mean-field and effective-medium approximations 
are also indicated by the dashed and solid lines, respectively. 

where z is the dynamical critical exponent. 

The numerical estimates of the autocorrelation time r*^ are plotted in Fig.^Jfor several 
dilutions. They show a growth of the autocorrelation time with the lattice size which 
becomes dramatic as the bond concentration increases and a behaviour compatible with 
a power law of the system size when p decreases, as expected since the dilution softens 
the transition and thus reduces the dynamical exponent z. Nevertheless, it is not possible 
to distinguish precisely the two regimes on a plot of the autocorrelation time versus the 
lattice size. Here, we may locate approximately the boundary between the two regimes 
around - slightly above p = 0.68. Indeed, the autocorrelation time at p = 0.68 is very 
well fitted with a power-law for all lattice sizes smaller than L = 30. Above, the data 
display a downward bending that can be explained by a correction to the power-law 
behaviour but not by an exponential prefactor (the bending would be upward). On the 
other hand, for p = 0.84 it is not possible to find any set of three consecutive points 
that could be fitted by a power-law: the autocorrelation time clearly grows faster than 
a power-law. Using two successive lattice sizes Li and L2 > Li, we defined an effective 
dynamical exponent 



which is expected to reach a finite value for continuous transitions and to diverge for 
first-order ones. The data, plotted in Fig. 1^01 again do not lead to any sound estimate of 
the location of the tricritical point. Nevertheless, the transition again definitely remains 
continuous up to the bond concentration p = 0.68. For higher concentrations, the data 
show an increase of the dynamical exponent with lattice size, but it is not possible to 
state unambiguously whether they develop a divergence or not. We also notice that the 




(32) 
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Figure 19: Autocorrelation time of the energy r*^ with respect to the lattice size at the (pseudo-) 
transition temperature. The curves correspond to different bond concentrations p 
(from bottom p — 0.28 to the top p — 1.00 in steps of 0.04). The results shown here 
all follow from MC simulations using the Swendsen-Wang algorithm. 

necessary finite number of iterations leads to an underestimate of and thus of z for 
bond concentrations close to p = 1 at large lattice sizes (this is particularly clear in 
Fig. I^njfor the size L = 13— 16). Multi-bondic simulations were thus needed in this case 
to improve the measurement of thermodynamic quantities when p is close to 1. 

Another approach is provided by the behaviour of the order-disorder interface tension. 
Numerically, the interface tension can be estimated from the probability distribution 
P{e) of the energy. One has 

oc e-^^^ = e-2-o.d.i^-\ (33) 

Indeed, the free-energy barrier can be related to the ratio of the (equally high) probabil- 
ities of the ordered and disordered phases (corresponding to the two peaks) and of the 
mixed phase regime involving two interfaces^ and which corresponds to the bottom of 
the gap between the two peaks. We started from the effective transition temperatures 
estimated from the maxima of the magnetic susceptibility. At this temperature, the sta- 
tistical weight of the ordered and disordered phases are comparable so the height of the 
peaks is very different. In order to define the interface tension, we reweighted the time 
series of the simulations to the (close) temperature for which the two peaks have equal 
heights. The order-disorder interface tension is plotted against the inverse of the lattice 
size at the transition temperature in the upper part of Fig. It shows undoubtedly a 
vanishing of the interface tension for p = 0.56, and presumably for p = 0.76 (not shown 

®Due to the employed periodic boundary conditions only an even number of interfaces can occur for 
topological reasons. 
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Figure 20: Effective dynamical exponent (SW algoritlim) with respect to tlic smaller lattice size 
at the transition temperature. The curves correspond to different bond concentra- 
tions p (from bottom p = 0.28 to the top p= 1.00 in steps of 0.04). 
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Figure 21: Probability distribution of the energy at the temperature for which the two peaks 
have equal heights. The two plots correspond to two different bond concentrations: 
p = 0.56 on the left (SW algorithm, increasing sizes L — 25, 30, 35, 40, 50, 64, 
and 96) and p — 0.84 on the right (multi-bondic simulations, sizes L = 16, 20, and 
25). The order-disorder interface tension do.d. = ln(Pmax/fmin)/(2i^~^) is plotted 
against in the upper part of the figure. 
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here) also, being a clear indication of a disorder induced second-order transition. On the 
other hand, for p = 0.84 the interface tension seems to converge towards a finite (but 
very small?) value in the thermodynamic limit, which can be taken as a signal for the 
persistence of the first-order nature of the transition in the pure case at p = 1 down to 
this dilution. 

As a consequence, we are led to the conclusion that the tricritical point is presumably 
located between p = 0.68 and p = 0.84, the upper bound corresponding to the observa- 
tion of an exponential growth of the autocorrelation time and the lower to a constant 
dynamical exponent and the vanishing of the latent heat (both values of p are indicated 
in the previous Figs. 1191 and l20() . However, one cannot unambiguously prove by numerical 
simulations on finite systems that what we identified as a second-order phase transition 
is not a weak first-order phase transition with a correlation length larger than L = 128, 
or that the fast growth of the autocorrelation time for p > 0.84 is not a cross-over to a 
power-law regime at larger system sizes. 

6 Critical behaviour 

6.1 Leading behaviour and critical exponents 

We now concentrate on the second-order regime only, i.e., on p < 0.68 where we per- 
formed an investigation of the universality class at the disorder fixed point. The critical 
exponents are computed using the finite-size scaling behaviour of the physical quantities 
(Eqs. (|llj) - (|14)) ) at the effective transition temperature (3c{L,p). In the usual renormal- 
isation group scheme for disordered systems, the renormalisation flow is subject to the 
influence of three flxed points describing respectively the pure system, the random sys- 
tem and the percolation transition. The scaling behaviour is thus expected to display 
large corrections resulting in a cross-over to a unique universal behaviour at large lattice 
sizes. According to this scheme, the exponents which are measured are expected to be 
(apparently) concentration dependent. In the previous sections (see, e.g., Fig. [T7|l. the 
corrections to scaling for the transition temperature have been observed to be weaker 
for the bond concentration p = 0.56. This behaviour is illustrated, e.g., in Fig. 1221 where 
the cross-over effect reflects in the bending of the curves Pma.x{L,p)J vs. for three 
dilutions p = 0.32, p = 0.56, and p = 0.80. The curve at p = 0.56, on the other hand, 
is almost flat. The corresponding data for the three main dilutions in the second order 
regime, p = 0.44, p = 0.56, and p = 0.68, are then plotted against L^^/'^ on the right 
part. Although the value of v is not yet known, we anticipate here the later result, using 
already the "to-be-determined-exponent". Again, the curve at p = 0.56 has an almost 
vanishing slope. As a consequence, we decided to make further large-scale Monte Carlo 
simulations at this concentration p = 0.56 up to the lattice size L = 96. To monitor the 
effects of the competing fixed points, we also made additional large-scale simulations for 
the concentrations p = 0.44 (towards the percolation transition) and p = 0.68 (towards 
the regime of first-order transitions) up to the lattice sizes L = 128 and L = 50, respec- 
tively (size limitations at these concentrations are linked to the discussion of Sect.OJ. 
In Fig. 1231 the finite-size scaling behaviour of the maximum susceptibility, Xmax) the 
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Figure 22: Evolution of tiie size-dependent (pseudo-)critical coupling with the inverse system 
size for relatively small sizes on the left plots. The same on the right plot for the 
three main dilutions, where the data are by anticipation fitted to a linear relation 
/3max(i,p) = /3c(p) + aL~^/'^ + . . ., where our estimate for v (« 0.75) will be discussed 
later. The slope coefficient is slightly positive for p = 0.44, slightly negative for 
p = 0.68 and virtually zero at p = 0.56, where the corrections-to-scaling (at least for 
this quantity) appear to be the smallest. 
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Figure 23: Finite-size scaling behaviour of the susceptibility, the magnetisation and of 



(3L^^d\nrn/df3 at /3max (the quantities have been shifted in the vertical direction 
for the sake of clarity). The behaviour at small lattice sizes is presumably governed 
by the percolation fixed point (shown as dashed lines and characterised by exponent 
ratios j/i^ ~ 2.05 and P/i' ~ 0.475). Above a crossover length scale a new (random) 
fixed point is reached (shown by the solid lines, with exponent ratios j/v ~ 1.535, 
l/iy ~ 1.34, and P/i^ ~ 0.73, discussed in detail below). 



magnetisation at /3max and the derivative of In m with respect to the inverse temperature 
evaluated at /9max are plotted versus the system size on a log-log scale. These curves 
should give access to the exponents 7/1^, P/i^, and l/i^, respectively. The three main 
dilutions are represented. One clearly observes a crossover between two regimes. For 
small lattice sizes, the system is strongly influenced by the proximity of a perturbing fixed 
point while a different, unique fixed point, is apparently reached at large sizes, as revealed 
by the slopes which are at first sight independent of the dilution when the linear extent 
of the lattice reaches values of about L > 30. The most probable susceptibility X50% is 
shown in Fig. |^ and can also lead to estimates for 7/2^. According to the discussion 
given in Sect. El we expect that the most probable susceptibility is better described than 
the average susceptibility, for which there exists a significant contribution of rare events, 
and these rare disorder realizations might be poorly scanned if a too small number 
of samples is considered. This difficulty might be circumvented through the study of 
what we defined as X50% Eq. (|28j) . In the presence of multifractality, the universal 
behaviour of X50% should differ from that of x- Since such a peculiar behaviour does 
not occur in the case of a global quantity [2Sj, like we expect compatible values of 
7/1^ as deduced from X50% X- Observing the data plotted in Fig. in fact, confirms 
our previous analysis. It seems that X50% is less influenced by the crossover effects than 
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Figure 24: Log-log plot of the average susceptibility (open symbols) and the typical susceptibility 
(filled symbols), as defined in Eq. H28|l by X50% for the three principal dilutions 
studied, indicating that the asymptotic scaling regime sets in earlier for the latter 
quantity. 

the average Xmax- order to support this statement, we will present the results of fits 
of the susceptibility in two different tables for the two regimes and for the three main 
dilutions: 

t> At small lattice sizes, the behaviour of x^ax ^^'^ "^/3max in all three cases com- 
patible with the percolation exponents (7/i^)perco — 2.05 and (/3/i^)perco — 0.475 
shown in Fig. I^Slby the dashed hues. This seems to be true (particularly in the 
case of the susceptibility) over a wider range of sizes for p = 0.44 than for p = 0.56 
or p = 0.68. This observation is compatible with a stronger influence of the perco- 
lation fixed point when p = 0.44, which is closer to the percolation threshold than 
the two other dilutions. Surprisingly, the assumption of a percolation influence 
is absolutely not confirmed^ by the behaviour at small sizes of the third quantity 
of interest, {dlnm/d(3)j3^^^. Due to the involved differentiation with respect 
to inverse temperature, the identification with percolation quantities becomes less 
obvious, but we do not have any explanation for this strange result. In Table 
we try to point out the influence of the percolation fixed point. This is achieved by 
power-law fits between a fixed minimum size Luiin — 4 up to an increasing maxi- 
mum size -Lmax below the value -L = 30 which apparently marks the modification 
in the behaviour of the physical quantities under interest. We first observe that 
f3/u starts from a value very close to the percolation value, and second, that X50% 
has always a lower exponent (i.e. more distinct from the percolation value). 

'^The expected percolation exponent would be 1/z/ ~ 1.124 while the slope at small sizes is larger than 
in the random regime where it takes a value close to l/u 1.35. 
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Table 3: Exponents deduced from the finite-size scaling behaviour of Xmax ^^"^ X50% in 

vicinity of the percolation fixed point (small sizes). Recall the percolation value 
(7/^)pGrco — 2.05 for comparison. 







P 


= 0.44 


P 


= 0.56 


P 


= 0.68 






"f/v deduced from 


^ jv deduced from 


^ jv deduced from 




^max 


Xmax 


X50% 


Xmax 


X50% 


Xmax 


X50% 


4 


8 


2.015 


1.902 


2.098 


1.916 


2.211 


1.915 


4 


13 


1.984 


1.866 


2.034 


1.818 


2.132 


1.720 


4 


20 


1.954 


1.833 


1.973 


1.748 


2.051 


1.579 


4 


30 


1.924 


1.808 


1.913 


1.691 


1.974 


1.500 



Table 4: Exponents deduced from the finite-size scaling behaviour of Xmax ^^'^ X50% in 

vicinity of the random fixed point (large sizes). The largest size taken into account in 
the fits is imax = 128 for p = 0.44, 96 for p = 0.56, and 50 for p = 0.68. 



p = 0.44 p = 0.56 p = 0.68 

^jv deduced from ^jv deduced from ^ jv deduced from 





Xmax 


X50% 


Xmax 


X50% 


Xmax 


X50% 


20 


1.724 


1.672 


1.571 


1.579 


1.541 


1.412 


25 


1.711 


1.664 


1.543 


1.587 


1.479 


1.471 


30 


1.706 


1.669 


1.518 


1.596 


1.438 


1.539 


35 






1.500 


1.581 


1.447 


1.645 


40 


1.703 


1.679 


1.502 


1.587 


1.464 


1.675 


50 


1.695 


1.657 


1.506 


1.593 






64 


1.680 


1.659 











At large sizes, for each quantity considered here, the curves corresponding to the 
three dilutions in Figs. [521 and evolve, after a crossover regime whose exact 
location depends on the value of p, towards a presumably unique power-law be- 
haviour which seems to remain stable then (solid lines in Fig. I23|). We thus believe 
that we have reached large enough sizes in order to get reliable estimates of the 
random fixed 'point exponents. This is only a visual impression, since in fact the 
effective exponents are still subject to significant variations, especially for the ex- 
treme dilutions p = 0.44 and p = 0.68. Effective exponents ^/v, P/v, and 1/v 
may be defined from power-law fits of Xmaxi "^/3max5 a'^'^ L~^dlnfn/d/3 between 
an increasing minimum size, L^i^, and a maximum one, Lmax- The value -Lmax is 
kept to the maximum available value L = 128, 96, and 50 for p = 0.44, 0.56, and 
0.68, respectively, and the results for the susceptibility are presented in Table ID 
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Table 5: Linear fits for XmnK' ^Anax' L^^d\nrfi/dp at /3max, leading to finite-size estimates 
of the combinations of critical exponents 7/1^, P/i' and l/v. These results correspond to 
the three main dilutions, and they are extracted from the finite-size scaling behaviour 
of the quantities at the temperature where the maximum of the average susceptibility 
is found by histogram reweighting. The results for dilutions p — 0.44 and p = 0.68 are 
less stable than for p = 0.56, reflecting the role of the crossover. 



p 




Lmax 




error 




error 


\lv 


error 


7/j/ + 2/3/j/ 


0.44 


30 


128 


1.706 


0.006 


0.544 


0.005 


1.395 


0.006 


2.794(16) 




40 




1.703 


0.008 


0.552 


0.007 


1.381 


0.008 


2.807(22) 




50 




1.695 


0.010 


0.540 


0.009 


1.358 


0.010 


2.775(28) 




64 




1.680 


0.016 


0.534 


0.014 


1.357 


0.016 


2.748(44) 


0.56 


30 


96 


1.518 


0.011 


0.588 


0.010 


1.389 


0.011 


2.694(31) 




35 




1.500 


0.014 


0.592 


0.012 


1.362 


0.013 


2.684(38) 




40 




1.502 


0.016 


0.608 


0.015 


1.353 


0.016 


2.718(46) 




50 




1.506 


0.026 


0.645 


0.024 


1.330 


0.025 


2.796(74) 


0.68 


25 


64 


1.479 


0.021 


0.343 


0.015 


1.505 


0.021 


2.165(51) 




30 




1.438 


0.031 


0.344 


0.022 


1.453 


0.030 


2.126(75) 




35 




1.447 


0.047 


0.342 


0.033 


1.437 


0.046 


2.13(11) 




40 




1.464 


0.075 


0.547 


0.051 


1.379 


0.075 


2.56(18) 



We see there that X50% is again better behaved (more stable) than the average 
susceptibility. 

Since we are mainly interested in the randomness fixed point, we now concentrate on 
fits at large system sizes. An exhaustive summary (i.e. for all three dilutions under 
interest) of the results of the fits performed at dilutions p = 0.44, p = 0.56, and p = 0.68 
is presented in Table 13 The corresponding effective exponents are also plotted against 
L" in Fig. ESj These results show that the data analysis is much more complicated than 
our previous preliminary determination of exponents in Table HI Again, the crossover 
between percolation and random fixed point behaviours is visible through the variation 
of effective exponents and the data present large corrections-to-scaling. 

A precise determination of the magnetic exponents is quite difficult. Indeed, as can 
be seen in Fig. the effective critical exponents {'y/i^)eS and {P/i')eS do not converge 
towards p-independent limits when Lmin — > -^max- The cross-over effects on the thermal 
quantities are much smaller. Indeed, the effective critical exponent VeS is converging to 
a roughly p-independent limit when L^nin -^max- We can give the following estimates 
for 7/1/ and l/v : 



p = 0A4: (7/i/)efr ~ 1.68(2), (l/i.)efr ~ 1.36(2) 
p = 0.56: (7/i/)efr ^ 1.51(3), (l/z^)efr ~ 1.33(3) 
p = Om: (7/i/)efr ~ 1.46(8), (l/i^)efr =^ 1-38(8) 



(34) 
(35) 
(36) 
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Figure 25: Effective critical exponents -f/i^ and as computed from a power-law fit between 
Liiiin and L^i^^, with imax fixed to the maximum available value L = 128, 96 and 
50 for p = 0.44, 0.56, and 0.68, respectively. They are plotted against L^^. The 
thin solid line shows the percolation values and the shadow stripe corresponds to our 
estimate for the random fixed point values. In the case of the dilution p = 0.56, the 
value of 4- 7/1/ is also shown. 

simply corresponding to the last line of Table 13 i.e., to the largest studied value of 
Lmin, for each dilution. The value of /3/z/ on the other hand is definitely not stable 
and more subject to the competing influence of fixed points. For p = 0.44 for example, 
the estimate of {l3/i')es is relatively stable against variations of L^[^, with fitted values 
below 0.5, close to the expected value for the percolation transition (0.475). This is a 
quantitative indication that the system is probably still subject to cross-over caused by 
the percolation fixed point. In the case of p = 0.68, the estimate of {P/i')ee is very 
small, then suddenly increasing for Lmin = 64. These remarks are consistent with the 
renormalisation scheme described above. In order to help us to decide between the 
different effective values measured at the three dilutions, we use the scaling relation 
^ /v + 2(3 /v = D = 3 which is almost satisfied for the bond concentration p = 0.56 
only (shown in Fig. I25|) when taking into account the lattice sizes L > 50. For the 
bond concentrations p = 0.44 and p = 0.68, this scaling relation is not satisfied for any 
of the accessible values. One is thus led to conclude that the critical regime has not 
yet been reached for these concentrations, in spite of our efforts to go up to very large 
sizes. Remember also that the corrections-to-scaling were found to be the smallest at 
p = 0.56, so the asymptotic regime in neighbouring dilutions should be more difficult to 
reach. Figure 123 thus suggests to rely only on the values measured at dilution p = 0.56, 
extrapolated to Lmin — > 00, as shown in Fig. 12^1 where a dashed stripe emphasises such 
an extrapolation of the effective exponents measured at the largest sizes. The values of 
(7/i/)ofr and (l/i^)efr are indeed stable in the regime L > 35. We may thus have reliable 
estimates of the asymptotic values for these exponents, and a reasonable estimate for 
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Figure 26: Effective critical exponents ^/v, (3/u, and for the dilution p = 0.56 obtained 
from fits between Lmin and imax = 96 and extrapolated to Z/min — > 00. In this limit, 
the scaling relation 7/1/ + 2(3/y = D ]s nicely satisfied. 

P/v, ratifying the scaling relation. 

Using this extrapolation procedure, our final estimates of the critical exponents of the 
disorder induced random fixed point of the three-dimensional bond-diluted 4-state Potts 
model are the following values: 



7/1/ = 1.535(30), 
P/v =0.732(24), 
l/v = 1.339(25), 



(37) 
(38) 
(39) 



resulting from a linear extrapolation of the data points for L^i^ = 25, 30, 35, 40, 50, 
and 64 at p = 0.56. Note that since the data are correlated, we have kept the error of 
the last point. 



6.2 Corrections to scaling 

For the 3D disordered Ising model it is well known that the correction-to-scaling close 
to the random fixed point are strong (with a corrections-to-scaling exponent around 
u = 0.4). Let us assume here also the existence of an irrelevant scaling field g with 
scaling dimension yg = —uj < 0. The scaling expression for the susceptibility 



(40) 



expanded at fic (on a finite system the susceptibility is always finite) around the fixed 
point value g = Q, leads to the standard expression 
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Figure 27: Plot of the deduced from linear fits of x^^xl-^) = TcL'^/''{l + b^L-'^) in the range 



25 < L < 96 for p = 0.56. The exponents are treated as fixed parameters and the 
amplitudes are free. The base plane gives the ranges of variation of the exponents: 
1.25 < ^jv < 1.75 and < w < 5. The absolute minimum is at 7/1^ = 1.49, 
w = 3.88, but the valley is extremely flat in the w-direction. A cutoff at = 50 has 
been introduced in order to improve clarity of the figure. 



order to investigate this question for the 3D 4-state Potts model, we tried to fit the 
physical quantities for p = 0.56 as 



and similar expressions for m^^^^^ in the range L > 25 where the leading term was 
already fitted in the previous section, and the subleading correction is due to the first 
irrelevant scaling field. 

Since four-parameter non-linear fits are not stable, we preferred linear fits where the 
exponents are taken as fixed parameters but the amplitudes are free. In Fig. 1271 we 
show a 3D plot of the cumulated square deviation of the least-square linear fit, ^ ^ as 
a function of 7/z/ and uj. There is a clear valley which confirms that is close to 
1.5, but the valley is so fiat in the w-direction that there is no clear minimum to give a 
reliable estimation of the corrections-to-scaling exponent. The same procedure for fijv 
is illustrated in the next figure (Fig. I28j) . Again, there is no way to get a compatible 
corrections-to-scaling exponent from the three fits, but the leading exponents are indeed 
close to l3/v = 0.71 (and 1/v = 1.35). Of course the minima of do not exactly coincide 
with the data presented in the table which should correspond to w ^ 00. 

7 Conclusion 

We studied the three-dimensional bond-diluted 4-state Potts model by large-scale Monte 
Carlo simulations. The pure system undergoes a strong first-order phase transition. The 




(41) 
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Figure 28: Plot of the deduced from linear fits of m (the exponent is thus negative) in the 
range 25 < L < 96 for p = 0.56. In the base plane, the range of variation of the 
exponents is — 1 < —jS/v < —0.5 and < < 5, and the minimum is at (i/v = 0.85, 
uj = 0.135. 



numerical estimates of the dynamical exponent z and of the interface tension give evi- 
dences for the existence of a disorder-induced tricritical point for bond dilutions between 
p = 0.68 and p = 0.84 below which the transition is softened to second order. Very strong 
crossover corrections are observed up to lattice size L < 30 — 40. The regime of the ran- 
dom fixed point is best observed for the bond concentration p = 0.56. From the values 
of the ratios of exponents measured at that concentration, 

= 1.535(30), (42) 
P/u = 0.732(24), (43) 
= 1.339(25), (44) 

the following estimates of the critical exponents are derived: 

7 = 1.146(44), (45) 
P = 0.547(28), (46) 
u = 0.747(14). (47) 

Let us mention that these exponents are in reasonably good agreement with recent star- 
graph high-temperature expansions ^01 of this model which give 7 = 1.00(3). The 
value of is eventually safe with respect to the bound i' > 2/D = 0.6666... of the 
stability of the random fixed point. In the random fixed point regime, we are unable 
to extract from the numerical data any reliable correction-to-scaling exponent (linked 
to the possible appearance of irrelevant scaling fields), even though it is clear that such 
corrections cannot be ignored. 
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In some sense, the outcome of this time-consuming work is disappointing, since we 
were not able to reach the asymptotic regime where exponents in the second-order regime 
of the phase diagram become dilution-independent, since the corrections to scaling are 
too strong, and since the tricritical point was not located with precision. We believe that 
this is due to the extreme difficulty of the problem and not to an unadapted approach. 
Perhaps we were too ambitious, but we have the feeling that the final values given for 
the critical exponents are reliable enough and should not be contradicted in the future 
by similar studies. 
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